function [SimResult,StationNum,Pos,CNR,InterferencePos]=Simulation(Monitor_Area, ...
    InterferenceNum, InterferenceSourcePower,InterferenceDistState, InterferencePos, InterferenceType,...
    StationDistState, StationInter,StationNum,StationPos,SimTime)


Monitor_Area_lon_max= max(Monitor_Area(:,1));     %提取监测区域经纬度范围
Monitor_Area_lon_min= min(Monitor_Area(:,1));
Monitor_Area_lat_max= max(Monitor_Area(:,2));
Monitor_Area_lat_min= min(Monitor_Area(:,2));
SimResult=1;

if StationDistState==1                %格网部署，StationInter有效
    Station_lon=Monitor_Area_lon_min:StationInter:Monitor_Area_lon_max;
    Station_lat=Monitor_Area_lat_min:StationInter:Monitor_Area_lat_max;
    StationNum=length(Station_lon)*length(Station_lat);
    Station_H=rand(StationNum,1);
    [StationPos_Lon,StationPos_Lat] = meshgrid(Station_lon,Station_lat);
    StationPos=[reshape(StationPos_Lon,[],1), reshape(StationPos_Lat,[],1), Station_H];

    in=inpolygon(StationPos_Lon,StationPos_Lat,Monitor_Area(:,1),Monitor_Area(:,2));
    StationPos_in=StationPos(in,:);
    StationNum=sum(sum(in));
    StationPos=StationPos_in;

    if StationInter <= 0
        SimResult=0;
    end

elseif StationDistState==2            %随机部署，StationNum有效
    StationNum_in = 0;
    StationPos = [];

    if StationNum <= 0
        SimResult=0;
    end

    while StationNum_in < StationNum
        StationPos_Lon=floor(rand(StationNum - StationNum_in,1) * (Monitor_Area_lon_max - Monitor_Area_lon_min) * 10) / 10 + ones(StationNum - StationNum_in,1)* Monitor_Area_lon_min;
        StationPos_Lat=floor(rand(StationNum - StationNum_in,1) * (Monitor_Area_lat_max - Monitor_Area_lat_min) * 10) / 10 + ones(StationNum - StationNum_in,1)* Monitor_Area_lat_min;
        StationPos_H=floor(rand(StationNum - StationNum_in,1) * 50)/10;
        in=inpolygon(StationPos_Lon,StationPos_Lat,Monitor_Area(:,1),Monitor_Area(:,2));
        StationPos=[StationPos;StationPos_Lon(in), StationPos_Lat(in), StationPos_H(in)];
        StationNum_in=StationNum_in+sum(in);
    end

elseif StationDistState==3            %位置导入，StationPos有效
    StationPos_Lon=StationPos(:,1);
    StationPos_Lat=StationPos(:,2);

    in=inpolygon(StationPos_Lon,StationPos_Lat,Monitor_Area(:,1),Monitor_Area(:,2));
    StationPos_in=StationPos(in,:);
    StationNum=sum(sum(in));
    StationPos=StationPos_in;

else
    SimResult=0;

end

%确定卫星的位置
GNSS_tx_num=4;
GNSS_tx_H=2e7;
SatellitePos=zeros(GNSS_tx_num,3);
for i=1:GNSS_tx_num                   %以监测区域经纬度范围随机生成
    SatellitePos(i,1)=rand(1,1)*(Monitor_Area_lon_max-Monitor_Area_lon_min)+Monitor_Area_lon_min;
    SatellitePos(i,2)=rand(1,1)*(Monitor_Area_lat_max-Monitor_Area_lat_min)+Monitor_Area_lat_min;
    SatellitePos(i,3)=GNSS_tx_H;
end

GNSS_tx_pos=SatellitePos;

% GPS信号发射，目前为GPS信号L1频点
global setting
setting.T_Sample = 0.02;      % 单位s
setting.freq_Sample = 5e6;    % 采样频率
setting.freq_IF = 2e6;        %信号中心频率（变换到中频）
setting.freq_CA = 1.023e6;    % CA码频率
setting.freq_D = 50;          % D码速率

T=setting.T_Sample;
%每个接收机接收到四个卫星所发射的信号SS=强度*伪码*数据*相位%
freq_Sample=setting.freq_Sample;
N=round(freq_Sample*T);
%1 卫星的接收到的强度，然后只考虑中途的衰减作为最后接收到的功率。
%卫星信号考虑大尺度衰落 Lbf=32.44+20lgF+20lgD  D的单位是km f的单位是MHz
f=1575.42e6;
Gain_dB=-func_txandre_fade(GNSS_tx_pos,StationPos,f)+14.3;    %卫星大尺度衰落
Gain_fd = sqrt(2*10.^(Gain_dB./10));
PRN_array_LL=[3,7,10,13];

%%
StationCNR_t=zeros(1,StationNum);
if InterferenceDistState==1                   %区域内随机设定
    InterNum_in = 0;
    InterPos = [];
    while InterNum_in < InterferenceNum
        InterPos_Lon=floor(rand(InterferenceNum - InterNum_in,1) * (Monitor_Area_lon_max - Monitor_Area_lon_min) * 10) / 10 + ones(InterferenceNum - InterNum_in,1)* Monitor_Area_lon_min;
        InterPos_Lat=floor(rand(InterferenceNum - InterNum_in,1) * (Monitor_Area_lat_max - Monitor_Area_lat_min) * 10) / 10 + ones(InterferenceNum - InterNum_in,1)* Monitor_Area_lat_min;
        InterPos_H=floor(rand(InterferenceNum - InterNum_in,1) * 50)/10;
        inter_in=inpolygon(InterPos_Lon,InterPos_Lat,Monitor_Area(:,1),Monitor_Area(:,2));
        InterPos=[InterPos;InterPos_Lon(inter_in), InterPos_Lat(inter_in), InterPos_H(inter_in)];
        InterNum_in=InterNum_in+sum(inter_in);
    end
    InterferencePos=InterPos;
end

CN0=zeros(25,StationNum);
for count=1:25

    %4 卫星接收到的多普勒频移
    for i=1:StationNum
        Doppler_range = -1000:20:1000;
        Doppler_index = randperm(length(Doppler_range),GNSS_tx_num);
        Doppler_array_LL(i,:)=Doppler_range(Doppler_index);%Doppler_array_LL = [0,0,0,0]
    end

    % 添加噪声 查一下实际过程中的信号和噪声的比值，然后往下做
    noise_dB=-204;        %-(mean(Gain_dB,2)+SNR);%噪声的功率
    for i=1:StationNum
        Power_noise= sqrt(2*10.^(noise_dB/10));     %高斯白噪声的功率。
        noise_re(i,:) = Power_noise*(randn(1,N)+1i*randn(1,N));
    end

    fs = 5e6;                 %干扰参数设置，目前为压制式单音干扰
    Ts=1/fs;
    point = (0:N-1)*Ts;
    t = 0:Ts:T-Ts;            %采样点
    f_STJ=2000000;            %干扰频率为卫星信号中心频率2e6
    jam=cell(InterferenceNum,1);
    for j=1:InterferenceNum
        [jam{j,1}]=-func_txandre_fade(InterferencePos(j,:),StationPos,f)+InterferenceSourcePower(j);     %干扰大尺度衰落
    end

    %5 做出基础信号
    temp_s = func_Gps_Signal_Generation_shift_001(PRN_array_LL);

    for m=1:StationNum
        SignalDatas = zeros(N,1);
        for n=1:4                  %4颗卫星
            fr = setting.freq_IF + Doppler_array_LL(m,n);    % 载波频率
            carr = cos(2*pi*fr*point);
            Signal=(Gain_fd(m,n).*temp_s(n,:).*carr).';
            SignalDatas = SignalDatas + Signal;
        end
        gpsSignal = SignalDatas;
        GpsSignal_H0 = gpsSignal +(noise_re(m,:)).';       %无干扰信号

        %7 添加干扰信号
        ud=rand(1,N)*2*pi;%uniformly distributed，产生0到2π的n个点的均匀分布
        position_of_STJrp=randperm(N,1);%产生1~n中1个随机数，不重复
        STJrp=ud(position_of_STJrp);%得到一个随机相位random phase，服从(0,2π)均匀分布
        GpsSignal_H1=GpsSignal_H0;
        t1=2e-6;

        for n=1:InterferenceNum
            jam_dB=[jam{n,1}];
            if InterferenceType(n)==1                      %单音干扰
                y_STJ1=sqrt(2*10^((jam_dB(m,1)+10*log10(0.8))/10))*cos(2*pi*f_STJ*t+STJrp);  %主径功率80%
                x1=y_STJ1.';
                y_STJ2=sqrt(2*10^((jam_dB(m,1)+10*log10(0.2))/10))*cos(2*pi*f_STJ*(t-t1)+STJrp); %第二径信号
                x2=y_STJ2.';
                GpsSignal_H1=GpsSignal_H1+x1+x2;      %信号+噪声+两径干扰

            elseif InterferenceType(n)==2                  %多音干扰
                f_MTJ=[4e5,5e5,6e5,7e5];
                position_of_MTJrps=randperm(N,4);%产生1~n中4个随机数，不重复，position_of_rps是1×4
                MTJrps=ud(position_of_MTJrps);%得到随机相位random phase，服从(0,2π)均匀分布，rps是1×4
                intermediate=rand(1,4);%生成4个0~1之间的随机数
                amp_per_comp1=intermediate*(2*10^((jam_dB(m,1)+10*log10(0.8))/10))/sum(intermediate);
                amp_per_comp2=intermediate*(2*10^((jam_dB(m,1)+10*log10(0.2))/10))/sum(intermediate);
                %把intermediate里面的每个数都变为原来的(2*10^(JNR/10))/sum(intermediate)，这样它们的和就是2*10^(JNR/10)
                y_MTJ=0;
                for ic=1:4
                    component1=sqrt(amp_per_comp1(ic))*cos(2*pi*f_MTJ(ic)*t+MTJrps(ic));
                    component2=sqrt(amp_per_comp2(ic))*cos(2*pi*f_MTJ(ic)*(t-t1)+MTJrps(ic));
                    y_MTJ=y_MTJ+component1.'+component2.';
                end
                GpsSignal_H1=GpsSignal_H1+y_MTJ;      %信号+噪声+两径干扰

            elseif InterferenceType(n)==3                  %线性扫频干扰
                B=2.046e6;         %扫频带宽
                T_unit=10e-5;      %扫频周期
                k=B/T_unit;        %线性扫频频率
                t_unit=0:Ts:(T_unit-Ts);
                num_unit=ceil(N/length(t_unit));%向上取整
                y_LFM1=[];
                y_LFM2=[];
                position_of_LFMrp=randperm(N,1);%产生1~n中1个随机数，不重复
                LFMrp=ud(position_of_LFMrp);%得到一个随机相位random phase，服从(0,2π)均匀分布
                for count_unit=1:num_unit
                    y_LFM_1period1 = sqrt(2*10^((jam_dB(m,1)+10*log10(0.8))/10))*cos(pi*k*t_unit.^2+LFMrp);
                    y_LFM_1period2 = sqrt(2*10^((jam_dB(m,1)+10*log10(0.2))/10))*cos(pi*k*(t_unit-t1).^2+LFMrp);
                    y_LFM1=[y_LFM1,y_LFM_1period1];
                    y_LFM2=[y_LFM2,y_LFM_1period2];
                end
                y_LFM1=y_LFM1(1:N);
                y_LFM2=y_LFM2(1:N);
                GpsSignal_H1=GpsSignal_H1+y_LFM1.'+y_LFM2.';   %信号+噪声+两径干扰

            elseif InterferenceType(n)==4                  %脉冲干扰
                tao=0.1;%占空比
                A_PJ1=sqrt(10^((jam_dB(m,1)+10*log10(0.8))/10)/tao);%PJ的幅度
                A_PJ2=sqrt(10^((jam_dB(m,1)+10*log10(0.2))/10)/tao);%PJ的幅度
                T_PJ=5e-5;%PJ的周期
                num_p=ceil(T/T_PJ);%脉冲个数
                n_1period=ceil(N/num_p);%一个周期内的采样点数
                y_PJ_1period1=[A_PJ1*ones(1,round(n_1period*tao)),zeros(1,round(n_1period*(1-tao)))];
                y_PJ_1period2=[A_PJ2*ones(1,round(n_1period*tao)),zeros(1,round(n_1period*(1-tao)))];
                y_PJ1=[];
                y_PJ2=[];
                for count_p=1:num_p
                    y_PJ1=[y_PJ1,y_PJ_1period1];
                    y_PJ2=[y_PJ2,y_PJ_1period2];
                end
                y_PJ1=y_PJ1(1:N);
                y_PJ2=y_PJ2(1:N);
                GpsSignal_H1=GpsSignal_H1+y_PJ1.'+y_PJ2.';

            end

        end

        CN0(count,m)=func_Zaizaobi(GpsSignal_H1,3)+12;    %有干扰载噪比

    end

end

for x=1:StationNum
    sum_CN=0;
    for y=1:count
        sum_CN=sum_CN+CN0(y,x);
    end
    StationCNR_t(1,x)=sum_CN/count;
    if StationCNR_t(1,x)<=10               %载噪比过低接收机失效
        StationCNR_t(1,x)=0;
    elseif isnan(StationCNR_t(1,x))
        StationCNR_t(1,x)=0;
    end
end

Pos=cell(SimTime,1);
CNR=cell(SimTime,1);
[Pos{1,1}]=StationPos;
[CNR{1,1}]=StationCNR_t;
if SimTime>1
    for z=1:SimTime-1
        position=zeros(StationNum,3);
        cn=zeros(1,StationNum);
        for a=1:StationNum
            if StationPos(a)~=0
                cn(a)=StationCNR_t(a)+(rand(1,1)-0.5);
            else
                cn(a)=StationCNR_t(a);
            end
            for u=1:3
                position(a,u)=StationPos(a,u)+(0.02*rand(1,1)-0.01);
            end
        end
        [Pos{z+1,1}]=position;
        [CNR{z+1,1}]=cn;
    end
end

if SimResult~=0
    SimResult=1;
end

end